Lecture 9

Spatio-temporal models

Sara Martino
sara.martino@ntnu.no

Dept. of Mathematical Science, NTNU

Janine Illian
Janine.Illian@glasgow.ac.uk

University of Glasgow

Jafet Belmont
Jafet.Belmont@glasgow.ac.uk

University of Glasgow

Motivation

  • Many real-world data sets vary in both space and time.

    • Examples: rainfall, temperature, air pollution, crop yield, disease spread.

  • We often want to model how similarity (correlation) changes across both dimensions.

The good news 👍

  • the framework works with spatio-temporal data as well

The bad news 👎

  • fully spatio-temporal models are complex

  • model fitting can take a long time

  • different types of spatio-temporal data structures lead to different types of complexities

Complexity of spatio-temporal models

  • additional dependencies: in space and time

  • spatial and temporal behaviour can be independent on each other, or dependent: i.e. properties of the spatial structure vary through time (or vice versa)

Several options are available in inlabru

  • Areal Model

  • Geostatistical and Point Process models

The Spatio-Temporal Framework

We define a stochastic process:

\[ Z(s,t), \quad s \in \mathcal{S} \subset \mathbb{R}^d, \quad t \in \mathcal{T} \]

The covariance function:

\[ C((s_1,t_1),(s_2,t_2)) = \text{Cov}[Z(s_1,t_1), Z(s_2,t_2)] \]

Our goal: specify \(C(\cdot)\) to capture spatial, temporal, and joint dependencies.

Separable and non-separable models

Separable models \[ C((s_1,t_1),(s_2,t_2)) = C_S(s_1,s_2) \times C_T(t_1,t_2) \]

Interpretation:

  • Space and time act independently.

  • Spatial correlation not dependent on time lag.

Advantages: 😀

  • Simpler estimation.

  • Fewer parameters \(\rightarrow\) Easier computation.

Disadvantages: 😔

  • Unrealistic for evolving or propagating phenomena.

Non-Separable models

Separable and non-separable models

Separable models \[ C((s_1,t_1),(s_2,t_2)) = C_S(s_1,s_2) \times C_T(t_1,t_2) \]

Interpretation:

  • Space and time act independently.

  • Spatial correlation not dependent on time lag.

Advantages: 😀

  • Simpler estimation.

  • Fewer parameters \(\rightarrow\) Easier computation.

Disadvantages: 😔

  • Unrealistic for evolving or propagating phenomena.

Non-Separable models

\[ C((s_1,t_1),(s_2,t_2)) = f(|s_2-s_1|,|t_2-t_1|) \] Advantages: 😀

  • More realistic for dynamic physical processes.

  • Can model propagation or decay that varies over time.

Disadvantages: 😔

  • More parameters → complex estimation.

  • Need more data for estimation

  • Interpretation may be less straightforward.

We focus on separable models but some non-separable ones are available in inlabru1

Areal Data

Spatio temporal models for areal data

\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988 (simulated data!)

Space

Spatio temporal models for areal data

\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988

Time

Spatio temporal models for areal data

\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988

Time

Spatio temporal models for areal data

The observation model \[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]

We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)

  1. Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.

Model 1 - Implementation

cmp_1 = ~ Intercept(1) +
  space(id_area, model = "besag", graph = g, scale.model = T) +
  time(year, model  = "rw2", scale.model = T)

lik = bru_obs(formula = Y~.,
              data = out,
              family = "poisson",
              E = E )
fit_1 = bru(cmp_1, lik)

Model 1 - Implementation

cmp_1 = ~ Intercept(1) +
  space(id_area, model = "besag", graph = g, scale.model = T) +
  time(year, model  = "rw2", scale.model = T)

lik = bru_obs(formula = Y~.,
              data = out,
              family = "poisson",
              E = E )
fit_1 = bru(cmp_1, lik)

Spatio temporal models for areal data

The observation model \[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]

We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)

  1. Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.

  2. Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)

    Here the spatial fields are different for each year, they are independent on each other.

Model 2 - Implementation

cmp_2 = ~ Intercept(1) + 
  space(id_area, model = "besag", graph = g, scale.model = T,
                               replicate = id_time) +
  time(year, model  = "rw2", scale.model = T)

fit_2 = bru(cmp_2, lik)

Model 2 - Implementation

cmp_2 = ~ Intercept(1) + 
  space(id_area, model = "besag", graph = g, scale.model = T,
                               replicate = id_time) +
  time(year, model  = "rw2", scale.model = T)

fit_2 = bru(cmp_2, lik)

Spatio temporal models for areal data

The observation model

\[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]

We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)

  1. Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.

  2. Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)

    Here the spatial fields are different for each year, they are independent on each other.

  3. Model 3: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \[ u_{st} = \phi u_{st-1} + w_{st}, \text{ with } w_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T \] Here we are adding more temporal dependence in the spatial term, using a AR1 type model Here the spatial fields are different for each year, they are independent on each other.

Model 3 - Implementation

cmp_3 = ~ Intercept(1) + 
  space(id_area, model = "besag", graph = g, scale.model = T,
                               group = id_time, control.group = list(model = "ar1")) +
  time(year, model  = "rw2", scale.model = T)

fit_3 = bru(cmp_3, lik)

Model 3 - Implementation

cmp_3 = ~ Intercept(1) + 
  space(id_area, model = "besag", graph = g, scale.model = T,
                               group = id_time, control.group = list(model = "ar1")) +
  time(year, model  = "rw2", scale.model = T)

fit_3 = bru(cmp_3, lik)

Spatio temporal models for areal data

The observation model \[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]

  1. Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)

  2. Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)

  3. Model 3: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \[ u_{st} = \phi u_{st-1} + w_{st}, \text{ with } w_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T \]

Compare

  • Differences between Model 2 and 3 are not evident when estimating the past, but they would give different predictions

  • Model 3 in this case seem to better fit the data (more about these measures later)

    model      DIC     WAIC
1 Model 1 15562.21 19898.87
2 Model 2 15562.21 19898.87
3 Model 3 15562.21 19898.87
                     mean    sd 0.025quant 0.975quant
Precision for space 12.48  1.23      10.21      15.06
GroupRho for space   0.90  0.01       0.87       0.92
Precision for time  47.61 21.87      18.06     102.17

Spatio temporal models for areal data in inlabru

  • Models can be stuck together using either

    • replicate - the different slices are independent but share the hyperparameter
    • group - different dependence structures are implemented

Types of group model

names(inla.models()$group)
[1] "exchangeable"    "exchangeablepos" "ar1"             "ar"             
[5] "rw1"             "rw2"             "besag"           "iid"


Note The group and replicate features can be used for more than space-time modeling!

Continuous space models

Continuous space models

(Simulated data)

Spatio temporal models for areal data

The observation model

\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]

We are going to see 3 different models for the linear predictor \(\eta_{st}\)

  1. Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.

Model 1 - Implementation

Define mesh and spde model

border_simplified = st_simplify(border, dTolerance = 25)

mesh = fm_mesh_2d(boundary = border_simplified,
  max.edge = c(30, 90), cutoff = 15,
   offset = c(-0.05, -0.25),
  crs = st_crs(border))

spde = inla.spde2.pcmatern(mesh = mesh, 
  prior.range = c(200, 0.5), # P(range < 100) = 0.5
  prior.sigma = c(1, 0.5)) # P(sigma > 1) = 0.5

Note

  • Coordinate system is in Km (this is recommended always!)

  • The size of the domain is ca \(630\times488 \text{ Km}^2\) we use this as a guideline to define the prior for the range.

Model 1 - Implementation

Fit the model

cmp = ~ Intercept(1) + 
  time(time, model = "rw2") +
  space(geometry, model = spde)

lik = bru_obs(formula =  Y~ .,
              data = dd)

fit_1 = bru(cmp, lik)

# plot time effect

p = fit_1$summary.random$time %>% 
  ggplot() + 
  geom_ribbon(aes(ID, 
                  ymin = `0.025quant`,
                  ymax = `0.975quant`), 
              alpha = 0.5) +
  geom_line(aes(ID, mean))

Model 1 - Implementation

Space-time Predictions

pxl = fm_pixels(mesh, mask = border)

pxl_time = fm_cprod(pxl, data.frame(time = seq(1,12)))

pred_1 = predict(fit_1, pxl_time, ~ data.frame(mu = Intercept + space + time,
                                      time_effect = time,
                                      space_effect = space))

Spatio temporal models for areal data

The observation model

\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]

We are going to see 3 different models for the linear predictor \(\eta_{st}\)

  1. Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.

  2. Model 2: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \(u_t(s)\sim\text{GF}(\rho,\sigma_u),\ t = 1,\dots,T\)

    Here spatial fields are different for each year, are independent and share range and sd.

Model 2 - Implementation

Fit the model

cmp = ~ Intercept(1) + 
  time(time, model = "rw2") +
  space(geometry, model = spde, 
        replicate = time)

lik = bru_obs(formula =  Y~ .,
              data = dd)

fit_2 = bru(cmp, lik)

# plot time effect

p = fit_2$summary.random$time %>% 
  ggplot() + 
  geom_ribbon(aes(ID, 
                  ymin = `0.025quant`,
                  ymax = `0.975quant`), 
              alpha = 0.5) +
  geom_line(aes(ID, mean))

Model 2 - Implementation

Space-time Predictions

pred_2 = predict(fit_2, pxl_time, ~ data.frame(mu = Intercept + space + time,
                                      time_effect = time,
                                      space_effect = space))

Spatio temporal models for areal data

The observation model

\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]

We are going to see 3 different models for the linear predictor \(\eta_{st}\)

  1. Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.

  2. Model 2: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \(u_t(s)\sim\text{GF}(\rho,\sigma_u),\ t = 1,\dots,T\)

    Here spatial fields are different for each year, are independent and share range and sd.

  3. Model 3: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \[ u_t(s) = \phi\ u_{t-1}(s) + \omega(s),\text{ with } \omega(s)\sim\text{GF}(\rho,\sigma_u) \]

    Here we are adding more temporal dependence in the spatial term, using a AR1 type model Here the spatial fields are different for each year, they are independent on each other.

Model 3 - Implementation

Fit the model

h.spec <- list(rho = list(prior = 'pc.cor1', 
                          param = c(0, 0.9)))
cmp = ~ Intercept(1) + 
  time(time, model = "rw2") +
  space(geometry, model = spde, 
        group = time,
        control.group = list(model = 'ar1', 
                             hyper = h.spec))

lik = bru_obs(formula =  Y~ .,
              data = dd)

fit_3 = bru(cmp, lik)

# plot time effect

p = fit_3$summary.random$time %>% 
  ggplot() + 
  geom_ribbon(aes(ID, 
                  ymin = `0.025quant`,
                  ymax = `0.975quant`), 
              alpha = 0.5) +
  geom_line(aes(ID, mean))

Model 3 - Implementation

Space-time Predictions

pred_3 = predict(fit_3, pxl_time, ~ data.frame(mu = Intercept + space + time,
                                      time_effect = time,
                                      space_effect = space))

Adding covariates to space-time models

Non-separable models in inlabru

Summary

  • Separable models are easily implemented using the replicate and group features